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Abstract 

The initial data for black hole collisions is constructed using a conformal-imaging approach and 
a new adaptive mesh refinement technique, a fully threaded tree (FTT). We developed a second- 
order accurate approach to the solution of the constraint equations on a non-uniformly refined high 
resolution Cartesian mesh including second-order accurate treatment of boundary conditions at the 
black hole throats. Results of test computations show convergence of the solution as the numerical 
resolution is increased. FTT-based mesh refinement reduces the required memory and computer 
time by several orders of magnitude compared to a uniform grid. This opens up the possibility of 
using Cartesian meshes for very high resolution simulations of black hole collisions. 



1 Introduction 

This paper deals with the construction of initial 
data for black hole collisions on a high resolution 
Cartesian adaptive mesh. The problem of black 
hole collisions is an important problem of the dy- 
namics of spacetime, and has applications to fu- 
ture observations of gravitational waves by grav- 
itational observatories on Earth and in space. 

The problem of black hole collisions is highly 
nonlinear and can only be solved numerically. A 
solution must be obtained within a large com- 
putational domain in order to follow the out- 
going gravitational waves far enough from the 
source. At the same time, very high resolution 
is required near the black holes to describe the 
nonlinear dynamics of spacetime. Integration of 
the collision problem on a three-dimensional uni- 
form mesh requires enormous computational re- 
sources, and this remains one of the major ob- 
stacles to obtaining an accurate solution. 



Adaptive mesh refinement (AMR) can be 
used to overcome these problems by introduc- 
ing high resolution only where and when it is 
required. AMR is widely used in many areas of 
computational physics and engineering. It has 
been applied in a more limited way in general 
relativity |l| . There are several types of AMR. 
In a grid-based AMR, a hierarchy of grids is cre- 
ated, with finer grids overlayed on coarser grids 
if a higher resolution is required || . An unstruc- 
tured mesh approach uses meshes consisting of 
cells of arbitrary shapes and various sizes ||. 
A cell-based approach to AMR uses rectangular 
meshes that are refined at the level of individual 
cells. This approach combines high accuracy of 
a regular mesh with flexibility of unstructured 
AMR f§. The new introduced fully threaded 
tree (FTT) structure, which we use here, leads 
to an efficient, massively parallel implementation 
of a cell-based AMR fl. 
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The first step in solving the black hole col- 
lision problem is to construct the initial data. 
The widely used conformal-imaging approach 
has been proposed in |f§ , Q , || and developed 
in P|, |lC||. Another approach to the construc- 
tion of initial data was recently proposed in 
jil| . Approaches for constructing initial data 
for certain specific cases of black hole collisions 
were proposed in ]i2| , |i3|| , |i4fl . The conformal- 
imaging approach .j|7[ . ]S[ . ]9|] , |lO| consists of 
constructing the extrinsic curvature (momentum 
constraints equations) using an imaging tech- 
nique and then solving a nonlinear elliptic equa- 
tion for the conformal factor (energy constraint) 
with an appropriate mirror-image boundary con- 
dition. This approach is adopted in this paper. 

A numerical technique for obtaining initial 
data for black hole collisions on a uniform Carte- 
sian grid using conformal-imaging approach is 
described in ||l0| . Two major problems with this 
approach mentioned in jl0| are the low resolu- 
tion of a uniform grid near black holes, and low- 
order accuracy and programming complexity of 
the inner boundary conditions at the black hole 
throats. In |iC|] , first-order accurate boundary 
conditions were implemented. The goal of this 
paper is to construct initial data for black hole 
collisions on a high resolution, Cartesian FTT 
adaptive mesh. In the process of realizing this 
goal, we found that the accuracy of the solution 
critically depend on the accuracy of the numeri- 
cal implementation of the inner boundary condi- 
tion. We developed a simple second-order accu- 
rate algorithm for boundary conditions to deal 
with this difficulty. 

The paper is organized as follows. The next 
Section || presents the formulation of the prob- 
lem and the equations solved. Section || de- 
scribes the FTT technique, the finite-difference 
discretization of the problem, and the numerical 
solution techniques. Section |J presents the re- 
sults of the solutions for various configurations 
of two black holes configurations and compares 
these with existing solutions. 

2 Formulation of the problem 

The ADM or 3+1 formulation of the equa- 
tions of general relativity works with the met- 
ric 7?„- and extrinsic curvature Kf- of three- 
dimensional spacelike hypersurfaces embedded 
in the four-dimensional space-time, where i, j = 
1, 2, 3, and the superscript ph denotes the phys- 



ical space. On the initial hypersurface, 7^ and 

Kfj must satisfy the constraint equations ||. 
The conformal-imaging approach assumes that 
the metric is conformally flat, 

7^0 4 7 y , (1) 

where 7^ is the metric of a background flat 
space. This conformal transformation induces 
the corresponding transformation of the extrin- 
sic curvature 

k',', o 2 . (2) 
With the additional assumption of 

trK = , (3) 
the energy and momentum constraints are 

V 2 <f>+^<f>- 7 K l3 KV =0 , (4) 

and 

DjK ij = , (5) 

respectively, where V 2 and Dj are the laplacian 
and covariant derivative in flat space. 

A solution to (|^) for two black holes with 
masses Mg, linear momenta Pg and angular mo- 
menta Sg, where 5 = 1, 2 is the black hole index, 
is § 

K ij {r) = Kf?{v)+K%%r) , (6) 

where 

2 / 1 

Kfp{v) =3^ —2 (P Stine j + Psjn Bti - 
8=1 

(7) 

- ng ti n St j) P St kn s )) 

and 

nl k \\ (8) 

ekjib s n s n 5 ,i)) . 

In (@) and (§, the comma in the subscripts 
separates the index of a black hole from the co- 
ordinate component indices and is not a symbol 
of differentiation, Rg = Mg/2 is the black hole 
throat radius, and 115 = (r — rg)/\r — rg\ is the 
unit vector directed from the center of the <5-th 
black hole r,5 to the point r. We work in units 
where G = 1, c = 1. 
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Fi glire 1' This figure shows schematically the computational domain used in the computations. The computational 
domain is a cubic box of size L. Two black holes with throat radii _Ri and R2 are positioned on the X axis in the XY 
plane (z = 0) at equal distances from the origin, O, of the coordinate system. 



The inversion-symmetric solution to (|^) can 
be obtained from (|^) by applying an infinite se- 
ries of mirror operators to (|^), as described in j9|. 
Note, that before applying the mirror operators 
to K^™ 3 , this term must be divided by 2 since 
the image operators will double its value. The 
series converges rapidly, and in practice only a 
few terms are taken. In this paper we take first 
five terms (for details see ]l5[). 

After the isometric solution for Kij is found, 
(^) must be solved subject to the isometry 
boundary condition at the black hole throats 



n\Di(j) = 



2R S 



(9) 



and the outer boundary condition <fi — > 1 when 
r — > 00. This boundary condition is represented 

by 1 



dr 



1 



(10) 



where r is the distance from the center of the 
computational domain to the boundary. 

3 Numerical method 

We work in the Cartesian coordinate system 
in the background flat space and find the so- 
lution within a cubic computational domain of 
size L. Figure [j] shows the schematic represen- 
tation of the computational domain. Centers of 
the throats of black holes are located in the XY 
plane z = on the line y = at equal distances 
from the origin. The space inside the black hole 
throats is cut out. The inner boundary condition 



(||) is imposed at the surface of the two spheres, 
and the outer boundary condition is imposed at 
the border of the computational domain. 

Figure || shows the cut of the computational 
domain through the z = plane and gives an 
example of an adaptive mesh used in computa- 
tions. The values of all variables are defined at 
cell centers. Figure || shows that coarse cells are 
used at large distances from the black holes, and 
the finest cells are used near the throats where 
the gradient of 4> is large. The grid is refined 
to achieve a desired accuracy of the solution as 
described below. 

There are three types of cells. The first type 
are internal cells that are actually used in com- 
putations (these cells are white in Figure ||) . The 
layer of boundary cells which is one cell wide 
along the outer border of the computational do- 
main is used to define the outer boundary con- 
ditions. The layers of boundary cells inside the 
throats, two cells wide, are used to define the 
inner boundary conditions. A cell is considered 
as located inside a throat if its center is located 
inside. Boundary cells are indicated as shaded 
on Figure [|. 

We solve (|j) as follows. Similar to we 
introduce a new unknown variable 



u = o — a 



where 



8=1 



lis 



r - r s \ 



(11) 



(12) 



The reason for using this transformation is that 
u varies slower than near the throats, and is 



Figure 2: A cut of a computational domain along the XY plane (z — 0) showing an adaptive mesh. The resolution 
increases near the black holes. Internal cells are unshaded while boundary cells are shown shaded. 



more c onve nient for numerical calculations (see 
Section 4.1). In Cartesian coordinates, (jij) then 
becomes 



V 2 u + F(u) = 



where 



and 



F(u) 



(l + cm) 



(13) 



(14) 



(15) 



Equation flU) is a nonlinear elliptic equation. 
Below we describe the numerical procedure of 
finding its solution. First consider a cell of size A 
which has six neighbors of the same size. Let us 
number this cell and its neighbors with integers 
from to 6, respectively. Then the discretized 
form of @ is 



ui + u 2 + u 3 + u A + it 5 + w 6 - 6i* 



+ F(u o ) = 
(16) 



A finite-difference form of ([l3]) is more compli- 
cated for cells that have neighbors of different 
sizes and may involve larger number of neigh- 
bors in order to maintain second order accuracy. 
This is described in Appendix |A|. In general, 
for every internal cell, the finite-difference dis- 
cretization may be written as 



f(u ,ui,u 2 , -,u„) = 



(17) 



where u±,...,u n are the values of u in n neigh- 
boring points chosen to represent the finite- 
difference stencil of a cell. In the particular case 



of a cell with all equal neighbors, / is defined by 
the left-hand side of (16|). 



We solve the set of (17) by the Newton 
Gauss-Seidel method |16[ ], that is, we obtain a 
new guess of Uq cw using Newton iteration with 
respect to the unknown i*o 



uo - /(no, ■•■) 



( dfjug,...) 

V du 



(18) 



Then we accelerate the convergence by using a 
successive overrelaxation (SOR) 



u% ew = uu^ ew + (1 - co)u 



(19) 



where uj is the overrelaxation parameter. For a 
simple case of all equal neighbors, ( |l8| ) can be 
written as 



new _ 

u o — u o 



(itl + 1i 2 + 1*3 + U 4 + 1*5 + 1*6 — 6li ) 

6 - A 2 ( ^hol) 

\ duo ) 

(A 2 F(u )) 



(20) 



6 - A 2 l ^P^- 

duo 



For stencils with non-equal neighbors the dis- 
cretization of equation ( p"3| ) is given in Appendix 
[A|, from which the expressions for the Newton 
Gauss-Seidel iterations in those cases can be 
written explicitly. 

We select the value of to from the interval 
[l,w ma J by the method described in |l7]]. The 
value of w max is initially set to uj max — 1.995. If 
the solution begins to diverge during the itera- 
tions, oj is reset to 1, tu max is decreased by 2%, 
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Figure 3: An example of the computation of outer boundary values. Outer boundary cells are the three cells located 
at the top of the figure. Cell 1 is an outer boundary cell, with center at point C. Point A is the point located one cell 
size away from the C on the line connecting point C with the center of the computational domain. Point A is located 
inside cell 2, which is an internal cell centered on point B. The value of 4> a t point A, to be used in the outer boundary 
condition, is interpolated from point B with second order accuracy using the derivatives evaluated at point B. These 
derivatives depend on 4> at Cell 2 and it's neighboring cells, including outer boundary cells. 



and the relaxation is continued allowing u> to in- 
crease up to the new u) max or until the solution 
starts to diverge again. Iterations are terminated 
when 



( p"l| ) with a 1 evaluated in point A. The finite- 
difference expression of (|l0]) can be written as 



U Q — Mo 
U 



< e 



(21) 



for all Mo, where e is a predefined small number. 
In this paper, we do not attempt to accelerate 
the iterations using, for example, a multigrid or 
other sophisticated techniques since the initial 
value problem must be solved only once. 

The procedure described above assumes that 
the values of u are known at all neighbors. For 
internal cells that are close to a boundary, these 
values are substituted with the values of u in 
boundary cells. Now we describe the procedure 
of assigning values of u to boundary cells. A sim- 
ilar technique was used for fluid flow simulations 
about complex bodies p8| . 

Figure || illustrates the process for the outer 
boundary. In this figure, cell 1 is a boundary 
cell. We need to define a new value u^f w in its 
center, point C. We find another point, A, which 
is located at a distance A (equal to the size of 
cell 1) from point C along the line that connects 
C with the center of the computational domain. 
The value of ua is found by the second order 
interpolation using the values of u in cell 2 and 
all of its neighbors. The interpolation involves 
old values of u in both internal cells and bound- 
ary cells. Then <\>a is computed using equation 



1 - 



M/2 



{r c - A/2) 



(22) 



and then solved for $q w . The value of Uq £W is 
then finally found using equation ([ll]) with or 1 
calculated at point C. 

Figure [| illustrates the process of defining u 
for the inner boundary. Cell 1 is located inside a 
throat, and we need to define a new value Uq ew 
in its center, point C. The point A is a point out- 
side the throat that has the same distance to the 
inner boundary as C, along the normal to the 
throat passing through point C. Let us denote 
the distance between C and A as Aac- Again, 
the value of ua is found by second-order inter- 
polation using old values of u in cell 2 and its 
neighbors, and 4>a is calculated using equation 
111). Then the boundary condition (O) becomes 



- 4> A _ <pA + 4> n c 



(23) 



which can be solved for <jf^ w ; R$ is the radius of 



the throat that contains point C. As before the 
value of Up™ is then finally found using equation 
(p"T|) with aT 1 evaluated at point C. 



After all inner and outer boundary points are 
defined, the next iteration (21) is performed for 
internal cells and so on, until the iterations con- 
verge. The advantage of the method described 
above is that the new method applies the same 
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hole throat 



Figure 4: An example of the computation of inner boundary values. Cell 1, centered on point C, is the inner boundary 
cell. Point A is the point located outside of the black hole throat on the line between the center of the black hole throat 
and point C that has the same distance to the boundary as point C. Point A is located inside cell 2, centered on point 
B. The value at point A, to be used in the inner boundary condition, is interpolated from point B with second order 
accuracy using the derivatives evaluated at point B. These derivatives depend on the value of Cell 2 and it's neighboring 
cells, including inner boundary cells. 



numerical algorithm to all cells, and that this al- 
gorithm is second-order accurate. The method 
described in jlO| required 77 different numerical 
stencils corresponding to different relative posi- 
tions of the boundary and interior cells, and was 
only first-order accurate. 

The computational domain used by the FTT 
is a cube of size L. It can be subdivided 
to a number of cubic cells of various sizes 
1/2, 1/4, 1/8,... of L. Cells are organized in 
a tree with the direction of thread pointers in- 
verted. These pointers are directed from chil- 
dren to neighbors of parent cells, as described 
in S . The most important property of an FTT 
data structure is that all operations on it, in- 
cluding tree refinement and derefinement, can 
be performed in parallel. A computer memory 
overhead of FTT is extremely small: two inte- 
gers per computational cell. All coding was done 
by using the FTTLIB software library [f) that 
contains functions for refinement, derefinement, 
finding neighbors, children, parents, coordinates 
of a cell, and performing parallel operations. 

We characterize an FTT mesh by the min- 
imum and maximum levels of leaves (unsplit 
cells) present in the tree, l m in and l ma x- We 
construct an adaptively refined mesh by start- 
ing with one computational cell representing the 
entire computational domain and then by subse- 
quently subdividing it by a factor of two until we 
reach the level l m i n . We find a coarse solution 
for two black holes at level i m , n using u = 1 as 
an initial guess. After this solution is obtained, 



we identify the regions that require more refined 
cells. These regions are then refined once, and 
a new converged solution is obtained. The old 
coarse solution is used as an initial guess for the 
new one. The procedure is repeated until the 
level of refinement reaches l max - 

4 Test results 

4-1 One Schwarzschild black hole: one- 
dimensional test 

Before presenting the results of three- 
dimensional test computations on an FTT mesh, 
we will discuss the accuracy of our numerical 
method using a simpler one-dimensional test 
problem. A single Schwarzschild black hole at 
rest has an analytic solution for the conformal 
factor 



r 



(24) 



where R is the throat radius and r is the dis- 
tance in the background space from the center 
of the black hole throat. The one-dimensional 
test was performed on a uniform grid in order to 
assess the influence of our treatment of bound- 
ary conditions ( p2| ) and (23), and the effects of 
changing the variable to a new variable u ( pd| ) 
on the accuracy of the solution. In these calcu- 
lations, we used a uniform one-dimensional grid 
consisting of n = 16, 32, and 64 cells. Cells 2 
through n.— l were interior cells. Cells 1 and n 



were the inner and outer boundary cells, respec- 
tively. The computations were performed for the 
grid size L — 10, throat radius R = 1, and us- 
ing convergence criterion e = 6 x 10~ 15 . The 
throat was located between the first and second 
cell centers at a distance Ar from the center of 
the border cell 1. Three values of Ar = 0.05A, 
0.5A, and 0.95A were considered, where A is the 
cell size. When Ar = 0.5A, the throat is located 
exactly in the middle between the points A and 
C in (^2|) , and the inner boundary condition ( |22| ) 
becomes second-order accurate regardless of the 
order of interpolation that is used for finding ua- 
For Ar = 0.05A and 0.95A, the overall accuracy 
of the inner boundary condition depends on the 
interpolation used for finding ua- 

Numerical solutions were obtained using 
three different methods: (a) solving the finite- 
difference form of (Q) for the original unknown 
variable (f> and using first order interpolation to 
find 4>a (the rest of the boundary condition pro- 
cedure was identical to that described in Sec- 
tion |^), (b) solving for the original variable <fi 
but using second order interpolation, and finally, 
(c) using both the second order interpolation 
and solving for the new unknown u as described 
m Section |. Table compares the numerical 
and analytical solutions for these three cases by 
showing the maximum relative deviation of the 
numerical solution from the analytical solution 
for the interior points of the grid. As can be seen 
from Table |, the accuracy varies with the grid 
resolution (n), the method of interpolation, and 
the choice of the unknown variable (4> or u). It 
also depends on the exact location of the throat 
relative to grid points (Ar/A). As we expect 
various relative locations of the throats relative 
to grid points in three-dimensional calculations, 
we need a numerical procedure that provides a 
second-order accuracy in all cases. 

Results using method (a) show that the ac- 
curacy of the solution using first order interpola- 
tion for 4>a is unacceptable. The third row of Ta- 
ble [l] shows that the accuracy does not improve 
with increasing number of cells. Results using 
method (b ) show that second order interpolation 
for 4>a (method (b )) leads to the overall second- 
order algorithm. The accuracy of the solution 
increases roughly by a factor of four when the 
grid resolution is doubled. Results for method 
(c) shows that the accuracy is further improved 
for small Ar. 

It is possible to give an analytical estimate 



for an error introduced by the numerical inner 
boundary condition ( |22| ) for the Schwarzschild 
black hole case. In this case, the general solu- 
tion of (0), limited at infinity, is 4> — c + b/r. 
Let us assume that the numerical outer bound- 
ary condition does not introduce any error, and 
the solution in interior points is found exactly. 
Thus, the numerical approximation of the inner 
boundary condition is a unique source of numer- 
ical error. Then the numerical solution would 
have the form <j> = 1 + b/r. The difference be- 
tween b and R in (53) then will determine the 
overall error in the solution. We can find b by 
substituting cj> = 1 + b/r into (p2|). The estimate 
of the relative error then is 



Relative error = 



R-b 
R 



Ar 
2R 



(25) 



The estimate of the error using ( p5| ) is given in 
Table |l| in brackets for the method (c). The 
comparison with the numerical error indicates 
that for method (c), the error in the solution is 
second-order and is determined by the accuracy 
of the inner boundary condition rather than by 
errors of numerical calculations for internal cells. 



4-2 Time- symmetric initial data for two 
black holes. 

Next, we consider the case of two black holes 
with Pg = Ss = that have masses Mi = 1, 
M2 = 2 and located (positions of the centers of 
their throats) at r x = (-4,0,0), r 2 = (4,0,0) 
with finite separation |ri — r2 1 =8. The size of 
the computational domain is L = 64. Numeri- 
cal solutions were obtained using FTT adaptive 
meshes with different increasing resolutions near 
the black hole throats. We characterize the reso- 
lution by specifying the minimum and maximum 
levels of cells in the tree, l m in an d Imax- The 
cell size at a given level I is A; = L ■ 2~ l . The 
computations were performed on meshes with 
lmin = 4,5,6 and l max = 6,7,8,9,10,11. The 
refinement criterion for this case was the require- 
ment that 




dy 



dz 



< 0.05 



1/2 



(26) 



in every cell, where partial derivatives in ( pq ) are 
determined by the numerical differentiation. We 
used e = 6 x 10~ 7 in (En) to terminate iterations. 
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Accuracy of one- dimensional computations. 



Method 


Ar/A 


n= 16 


n=32 


n=64 




0.95 


-2.9 x 10" 1 


-7.2 x 10" 2 


-1.8 x 10" 2 


(a) 


0.5 


-6.1 x 10" 2 


-1.7 x 10~ 2 


-5.0 x 10~ 3 




0.05 


-3.8 x lO" 1 


-4.6 x 10" 1 


-4.7 x 10" 1 




0.95 


: t : 

-2.9 x 10" 1 


n 

-7.0 x 10~ 2 


n 

-1.7 x 10~ 2 


(b) 


0.5 


-6.1 x 10~ 2 


-1.7 x 10~ 2 


-5.0 x 10~ 3 




0.05 


-3.0 x 10- 1 


-1.2 x 10" 1 


-4.2 x 10~ 2 




0.95 


-3.4(3.5) x 10" 1 


-8.7(8.8) x 10" 2 


-2.2(2.2) x 10" 2 


(c) 


0.5 


-7.4(9.8) x lO- 2 


-2.1(2.4) x 10" 2 


-5.7(6.1) x 10~ 3 




0.05 


-9.9(9.8) x 10" 4 


-1.4(2.4) x 10" 4 


-5.7(6.1) x 10" 5 



Table 1: Relative accuracy of the numerical solutions of (Q) for a Schwarzschild black hole obtained using three differ- 
ent methods (a,b,c), different resolutions (n), and different location of the throat relative to a grid (Ar). Numbers in 
parentheses for the method (c) is the accuracy estimate based on the consideration of the error introduced by the inner 
boundary condition (see Section ^l]) . 



The mirror-image symmetric analytic solu- 
tion for two time-symmetric black holes is given 
in Appendix [b| (for details of derivation see |L5[ ). 
Table ^ gives the comparison of the numeri- 
cal solutions with fixed l m in — 5 and varying 
Imax = 5 — 11 with the analytical solution. The 
table shows the maximum deviation of a numer- 
ical solution from the analytical one. It also 
shows the level of a cell where the maximum 
error was found. From the table we see that 
the accuracy of the solution increases approxi- 
mately linearly with increasing l max , and that 
the maximum error is located at maximum level 
of refinement near the throats. When we com- 
pare solutions obtained on different meshes on 
which the resolution was increased on all levels 
simultaneously (Table |^), we observe better than 
linear convergence, as it should be expected. 

The computations performed on an adaptive 
mesh allow us to save a significant amount of 
computational resources. For example, our solu- 
tion obtained on the l m in = 5, lmax = 11 adap- 
tive mesh used 6 x 10 5 computational cells. An 
equivalent uniform-grid computation with the 
same resolution near the throats would have re- 
quired using a 2048 3 uniform Cartesian grid with 
~ 8 x 10 9 cells. That is, in this case the compu- 
tational gain was ~ 10 4 . 

4-3 Two black holes with linear and angu- 
lar momenta 

Cook et al. JlCj considered the initial conditions 
for two black holes of equal mass Mi = M 2 = 2 



with non-zero linear and angular momenta (case 
A1B8 in Table 3 of @). In our coordinate 
system (Figure |l|), the components of linear 
and angular momenta of the black holes for the 
case A18B are Pi = -P 2 = (0,0,-14), and 
S x = (280, 280, 0) and S 2 = (0, 280, 280), respec- 
tively. The throats are located at ri = (4, 0, 0) 
and V2 — (—4, 0, 0) with a relative separation 
equal eight. We computed the A1B8 case us- 
ing the same size of the computational domain, 
L = 28.8 as that used in |10| , and using a series 
of refined meshes with increasing resolution near 
the throats, l m in = 5, I max = 5 — 11 . We used 
the same value of e as in Section |4^ but used a 
modified mesh refinement criterion 




Table ^ compares coarse solutions l m in = 5, 
lmax = 5 — 10 with the finest solution obtained 
on the l m i n — 5, l max = H mesh. It shows 
that both the maximum and the average devia- 
tion of the solutions decreases by more than two 
orders of magnitude when the maximum reso- 
lution near the throats is increased by a factor 
of 32. The solutions in [0 did not show im- 
provement with increasing resolution (see their 
Table |). Figure 5 shows the comparison of <f> 
on the line passing through the centers of the 
throats computed in this paper with resolution 
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Accuracy of two black hole time- symmetric computations. 



l>min 


Imax 


Max. error 


l err 


5 


- 5 


3.3 x 1(T 2 


5 


5 


- 6 


4.8 x 1(T 2 


6 


5 


- 7 


2.9 x 1(T 2 


7 


5 


- 8 


1.0 x 10~ 2 


8 


5 


- 9 


1.2 x 10~ 2 


9 


5 


- 10 


3.9 x 10~ 3 


10 


5 


- 11 


8.9 x 10~ 4 
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Table 2: The table shows the maximum relative deviation of the numerical solutions for the two Schwarzschild black 
holes (Section 4.2), and the level of cells l err where the maximum error is located for computations with different 
maximum resolution. 



Accuracy of two black hole time-symmetric computations. 



Imin 


Imax 


Max. error 


4 - 


8 


3.6 x 10~ 2 


5 - 


9 


1.2 x 10~ 2 


6 - 


10 


3.7 x 10~ 3 



Table 3: The table shows the maximum relative deviation of the numerical solutions for the two Schwarzschild black 



holes (Section 4.2) for computations where both minimum and maximum resolutions were increased simultaneously. 



Accuracy of computations of two black hole with non-zero linear and angular momenta. 



l>min 


Imax 


Max. error 


Avg. error 


5 


- 5 


1.1 x lO" 1 


1.1 x 10~ 2 


5 


- 6 


6.5 x 10~ 2 


9.0 x 10~ 3 


5 


- 7 


1.6 x 10~ 2 


1.8 x 10~ 3 


5 


- 8 


1.7 x 10~ 2 


2.1 x 10~ 3 


5 


- 9 


5.3 x 10~ 3 


6.4 x 10~ 4 


5 


- 10 


7.5 x 10~ 4 


6.7 x 10~ 5 



Table 4: Comparison of coarse solutions to the finest solution on the Imin = 5, Imax 
maximum relative deviation (Max. error) and the average relative deviation (Avg. e 
computed in Section 4.3. 



11 mesh. The table gives the 
or) of the numerical solutions 
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Figure 5: A comparison of <f> for the A1B8 case presented in JlfJ (their Figure 4.) (b) and our computations with 
resolution l m in = 5 and lmax = 8 (a). Differences are not visible on the scale of these graphs, except for the fact that we 
have better resolution near the holes which is why we have higher values close to the holes. 



Imin = 5, lmax = 8 with the results presented in 
JTof (their Figure 4). In our computations, finer 
cells cluster near the throats where the gradient 
in the solution is larger, whereas in [ jioj] , cells 
have the same size and are uniformly distributed 
in space. This combined with the overall second- 
order accuracy of our method is the reason why 
the adaptive mesh refinement solution improves 
when the mesh is refined (see Table |]) . 

5 Conclusions. 

In this paper, we applied a new adaptive 
mesh refinement technique, a fully threaded tree 
(FTT), for the construction of initial data for 
the problem of the collision of two black holes. 
FTT allows mesh to be refined on the level of 
individual cells and leads to efficient computa- 
tional algorithms. Adaptive mesh refinement is 
very important to the problem of black hole col- 
lisions because a very high resolution is required 
for obtaining an accurate solution. 

We have developed a second-order approach 
to representing both the inner boundary condi- 
tions at the throats of black holes and the outer 
boundary conditions. This allowed us to imple- 
ment an approach to the solution of the energy 
constraint that is formally second-order accu- 
rate. We presented results of tests for two black 
holes that demonstrated a good improvement of 
the accuracy of the solution when the numerical 
resolution was increased. 

The FTT-based AMR approach gives a gain 
of several orders of magnitude in savings of 
both memory and computer time (Section 



and opens up the possibility of using Cartesian 
meshes for very high resolution simulations of 
black hole collisions. A second-order boundary 
condition technique similar to that developed in 
this paper can be applied for the integration of 
initial conditions in time. We plan to use these 
techniques for time integration of the black hole 
collision problem. 
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A Appendix: Finite-difference for- 
mulas on the Fully Threaded 
Tree 



On the FTT mesh, we use the four different types 
of stencils shown in Figure ^| when a cell has zero, 
one, two, or three neighbors that are two times 
larger. These stencils involve nine neighbors of 
a cell and the cell itself. Let us introduce the 
vector of partial derivatives of u at the center of 
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c) 



d) 



Figure 6: Examples of the different stencils used for computing derivatives are shown for the 4 different cases: a) 3 
large neighbors, b) 2 large neighbors, c) 1 large neighbor and d) large neighbors 



cell 



y- and z-directions the weights are: 



Du k = (u 



du du du d 2 u d 2 u d 2 i 



dx' dy' dz' dx 2 ' dxdy ' dxdz ' 



dy 2 ' dydz ' dz 2 



(28) 



which includes the value of the function itself as 
the zeroth component. We can express the val- 
ues of Ui in i = 1 , . . . , 9 neighboring cells with 
the second-order accurate Taylor expansion us- 
ing Duk- This leads to a linear system of ten 
equations for the ten Duu unknowns. From this 
system, Duk can be expressed as a weighted sum 
of the values of Ui in the cell and its neighbors 



Du k = ^w l k i 



(29) 



7=0 



Due to the limited number of stencils encoun- 
tered in the FTT structure, this can be done 
once and for all and the weights can be stored in 
an array. 

For three bigger neighbors in the positive x-, 



wi = ^(-48,-14,54,-1,15, 



w 3 



W7 



1 



1,15,-11,-11,2), 



w 2 = —(-48,-1,15,-14,54, 



52A 



52A 



26A ; 



-1,15,-11,2,-11), 
(-48,-1,15,-1,15, 

-14,54,2,-11,-11), 
-(-4,14,-2,1,-15, 

1,-15,11,11,-2), 



™5 = ^(1,0,-1,0,-1,0,0,1,0,0), 
w 6 = -^(1,0,-1,0,0,0,-1,0,1,0), 



26A ; 



1 



■(-4,1,-15,14,-2, 
1,-15,11,-2,11), 



w s = —(1,0,0,0,-1,0,-1,0,0,1), 

Wg 



A 2 ' 
1 



26A 2 



(-4,1,-15,1,-15, 



(30) 



14,-2,-2,11,11). 
For two bigger neighbors in the positive x- and 
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y-directions the weights are: the other three cases, while all other derivatives 

j are given by the standard formulas for central 



Wi = (—48, —15, 57, —1, 15, differences on a uniform grid. Equation ( [L3| 

then can be expressed as 

-2,14,-12,-11,3), 

1 , Du 4 + Du 7 + Du 9 + F(u ) = . 

™2 = ^(-48,-1,15,-15,57, (33) 

-2,14,-12,3,-11), 

W3 = _L (0; o, 0, 0, 0, -1, 1, 0, 0, 0), B Appendix: Conformal factor of 



Wi= ("8. 15, -1,1, -15, 



TWO TIME-SYMMETRIC BLACK HOLES. 



A solution for the conformal factor of two timc- 



2, 14, 12, 11, 3), symmetric black holes (jT9), |^ can be written in 

1 ,.. _ , „ ■> n n i n n\ (31) the Cartesian coordinates as [ 15(1 

W5 = ^(! 5 0, -1,0, -1,0,0, 1,0,0), ' ^ 

- oo 

w 6 = ^a(l, 0,-1, 0,0, 0,-1, 0,1,0), cf>(r) = 1 + £ (F? + , (34) 
107 = 7— ^(-8, 1, -15, 15, -1, 



n=l 



28A 2V ' ' ' ' ' with 
2,-14,12,-3,11), 

1 



F?- 1 (iJi/Pif 1 ) for n odd; 



w 8 = -^{l, 0,0, 0,-1, 0,-1, 0,0,1), F™ = ^ F^ 1 (iWiV 1 ) , for n even; 

1 1 , for n = 0, 

W9 = —(-2,0,0,0,0,1,1,0,0,0). ) 

f F^ 1 (ifc/pSf 1 ) , for n odd; (35) 

For one big neighbor in the positive x direction F 2 " = ) p^ 1 {Ri/p 7 ^ 1 ) , for n even; 

the weights are: " ^ ' for „ = 0) 

»i = 3^ (-24, -8, 30, -1, 7, -1, 7, -6, -6, 2), where 

^3 = ^(0,0, 0,-1, 1,0, 0,0, 0,0), p'Sa 1 = K~ 1 -rs\, (36) 

w 3 = —(0,0,0,0,0,-1,1,0,0,0), a = 1,2, 5 = 1,2, r 5 are the positions of the 

centers of black hole throats, R$ are the throat 



radii, 



w 4 = ^(-6,8,0,1,-7,1,-7,6,6,-2), 

W5 = i(l, 0,-1, 0,-1, 0,0, 1,0,0), fJi(xr 1 )= for n odd; 

^ (32) xi l = I J 2 (x^ 1 ) , for n even; 

W6 = —(1,0,-1,0,0,0,-1,0,1,0), [ r , forn = 0, 

^ = -^(-2,0,0,1,1,0,0,0,0,0), fJslxT 1 ), fornodd; (37) 

^ x 2 l = ^ J^x'^ 1 ) , for n even; 

™8 = (1,0,0, 0,-1, 0,-1, 0,0,1), (r, for n = 0, 

w 9 = ^ (-2, 0,0, 0,0, 1,1, 0,0,0). and 

For zero big neighbors the mixed second order J<s(x) = ( — — ^) (x — r^) + rg 

derivatives are given by the same weights as for V l x — r<5 l / (38) 
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